/*
 * Copyright (c) 1999, 2013, Oracle and/or its affiliates. All rights reserved.
 * ORACLE PROPRIETARY/CONFIDENTIAL. Use is subject to license terms.
 *
 *
 *
 *
 *
 *
 *
 *
 *
 *
 *
 *
 *
 *
 *
 *
 *
 *
 *
 *
 */

package java.math;

/**
 * A class used to represent multiprecision integers that makes efficient
 * use of allocated space by allowing a number to occupy only part of
 * an array so that the arrays do not have to be reallocated as often.
 * When performing an operation with many iterations the array used to
 * hold a number is only reallocated when necessary and does not have to
 * be the same size as the number it represents. A mutable number allows
 * calculations to occur on the same number without having to create
 * a new number for every step of the calculation as occurs with
 * BigIntegers.
 *
 * @author Michael McCloskey
 * @author Timothy Buktu
 * @see BigInteger
 * @since 1.3
 */

import static java.math.BigDecimal.INFLATED;
import static java.math.BigInteger.LONG_MASK;

import java.util.Arrays;

class MutableBigInteger {

  /**
   * Holds the magnitude of this MutableBigInteger in big endian order.
   * The magnitude may start at an offset into the value array, and it may
   * end before the length of the value array.
   */
  int[] value;

  /**
   * The number of ints of the value array that are currently used
   * to hold the magnitude of this MutableBigInteger. The magnitude starts
   * at an offset and offset + intLen may be less than value.length.
   */
  int intLen;

  /**
   * The offset into the value array where the magnitude of this
   * MutableBigInteger begins.
   */
  int offset = 0;

  // Constants
  /**
   * MutableBigInteger with one element value array with the value 1. Used by
   * BigDecimal divideAndRound to increment the quotient. Use this constant
   * only when the method is not going to modify this object.
   */
  static final MutableBigInteger ONE = new MutableBigInteger(1);

  /**
   * The minimum {@code intLen} for cancelling powers of two before
   * dividing.
   * If the number of ints is less than this threshold,
   * {@code divideKnuth} does not eliminate common powers of two from
   * the dividend and divisor.
   */
  static final int KNUTH_POW2_THRESH_LEN = 6;

  /**
   * The minimum number of trailing zero ints for cancelling powers of two
   * before dividing.
   * If the dividend and divisor don't share at least this many zero ints
   * at the end, {@code divideKnuth} does not eliminate common powers
   * of two from the dividend and divisor.
   */
  static final int KNUTH_POW2_THRESH_ZEROS = 3;

  // Constructors

  /**
   * The default constructor. An empty MutableBigInteger is created with
   * a one word capacity.
   */
  MutableBigInteger() {
    value = new int[1];
    intLen = 0;
  }

  /**
   * Construct a new MutableBigInteger with a magnitude specified by
   * the int val.
   */
  MutableBigInteger(int val) {
    value = new int[1];
    intLen = 1;
    value[0] = val;
  }

  /**
   * Construct a new MutableBigInteger with the specified value array
   * up to the length of the array supplied.
   */
  MutableBigInteger(int[] val) {
    value = val;
    intLen = val.length;
  }

  /**
   * Construct a new MutableBigInteger with a magnitude equal to the
   * specified BigInteger.
   */
  MutableBigInteger(BigInteger b) {
    intLen = b.mag.length;
    value = Arrays.copyOf(b.mag, intLen);
  }

  /**
   * Construct a new MutableBigInteger with a magnitude equal to the
   * specified MutableBigInteger.
   */
  MutableBigInteger(MutableBigInteger val) {
    intLen = val.intLen;
    value = Arrays.copyOfRange(val.value, val.offset, val.offset + intLen);
  }

  /**
   * Makes this number an {@code n}-int number all of whose bits are ones.
   * Used by Burnikel-Ziegler division.
   * @param n number of ints in the {@code value} array
   * @return a number equal to {@code ((1<<(32*n)))-1}
   */
  private void ones(int n) {
    if (n > value.length) {
      value = new int[n];
    }
    Arrays.fill(value, -1);
    offset = 0;
    intLen = n;
  }

  /**
   * Internal helper method to return the magnitude array. The caller is not
   * supposed to modify the returned array.
   */
  private int[] getMagnitudeArray() {
    if (offset > 0 || value.length != intLen) {
      return Arrays.copyOfRange(value, offset, offset + intLen);
    }
    return value;
  }

  /**
   * Convert this MutableBigInteger to a long value. The caller has to make
   * sure this MutableBigInteger can be fit into long.
   */
  private long toLong() {
    assert (intLen <= 2) : "this MutableBigInteger exceeds the range of long";
    if (intLen == 0) {
      return 0;
    }
    long d = value[offset] & LONG_MASK;
    return (intLen == 2) ? d << 32 | (value[offset + 1] & LONG_MASK) : d;
  }

  /**
   * Convert this MutableBigInteger to a BigInteger object.
   */
  BigInteger toBigInteger(int sign) {
    if (intLen == 0 || sign == 0) {
      return BigInteger.ZERO;
    }
    return new BigInteger(getMagnitudeArray(), sign);
  }

  /**
   * Converts this number to a nonnegative {@code BigInteger}.
   */
  BigInteger toBigInteger() {
    normalize();
    return toBigInteger(isZero() ? 0 : 1);
  }

  /**
   * Convert this MutableBigInteger to BigDecimal object with the specified sign
   * and scale.
   */
  BigDecimal toBigDecimal(int sign, int scale) {
    if (intLen == 0 || sign == 0) {
      return BigDecimal.zeroValueOf(scale);
    }
    int[] mag = getMagnitudeArray();
    int len = mag.length;
    int d = mag[0];
    // If this MutableBigInteger can't be fit into long, we need to
    // make a BigInteger object for the resultant BigDecimal object.
    if (len > 2 || (d < 0 && len == 2)) {
      return new BigDecimal(new BigInteger(mag, sign), INFLATED, scale, 0);
    }
    long v = (len == 2) ?
        ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) :
        d & LONG_MASK;
    return BigDecimal.valueOf(sign == -1 ? -v : v, scale);
  }

  /**
   * This is for internal use in converting from a MutableBigInteger
   * object into a long value given a specified sign.
   * returns INFLATED if value is not fit into long
   */
  long toCompactValue(int sign) {
    if (intLen == 0 || sign == 0) {
      return 0L;
    }
    int[] mag = getMagnitudeArray();
    int len = mag.length;
    int d = mag[0];
    // If this MutableBigInteger can not be fitted into long, we need to
    // make a BigInteger object for the resultant BigDecimal object.
    if (len > 2 || (d < 0 && len == 2)) {
      return INFLATED;
    }
    long v = (len == 2) ?
        ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) :
        d & LONG_MASK;
    return sign == -1 ? -v : v;
  }

  /**
   * Clear out a MutableBigInteger for reuse.
   */
  void clear() {
    offset = intLen = 0;
    for (int index = 0, n = value.length; index < n; index++) {
      value[index] = 0;
    }
  }

  /**
   * Set a MutableBigInteger to zero, removing its offset.
   */
  void reset() {
    offset = intLen = 0;
  }

  /**
   * Compare the magnitude of two MutableBigIntegers. Returns -1, 0 or 1
   * as this MutableBigInteger is numerically less than, equal to, or
   * greater than <tt>b</tt>.
   */
  final int compare(MutableBigInteger b) {
    int blen = b.intLen;
    if (intLen < blen) {
      return -1;
    }
    if (intLen > blen) {
      return 1;
    }

    // Add Integer.MIN_VALUE to make the comparison act as unsigned integer
    // comparison.
    int[] bval = b.value;
    for (int i = offset, j = b.offset; i < intLen + offset; i++, j++) {
      int b1 = value[i] + 0x80000000;
      int b2 = bval[j] + 0x80000000;
      if (b1 < b2) {
        return -1;
      }
      if (b1 > b2) {
        return 1;
      }
    }
    return 0;
  }

  /**
   * Returns a value equal to what {@code b.leftShift(32*ints); return compare(b);}
   * would return, but doesn't change the value of {@code b}.
   */
  private int compareShifted(MutableBigInteger b, int ints) {
    int blen = b.intLen;
    int alen = intLen - ints;
    if (alen < blen) {
      return -1;
    }
    if (alen > blen) {
      return 1;
    }

    // Add Integer.MIN_VALUE to make the comparison act as unsigned integer
    // comparison.
    int[] bval = b.value;
    for (int i = offset, j = b.offset; i < alen + offset; i++, j++) {
      int b1 = value[i] + 0x80000000;
      int b2 = bval[j] + 0x80000000;
      if (b1 < b2) {
        return -1;
      }
      if (b1 > b2) {
        return 1;
      }
    }
    return 0;
  }

  /**
   * Compare this against half of a MutableBigInteger object (Needed for
   * remainder tests).
   * Assumes no leading unnecessary zeros, which holds for results
   * from divide().
   */
  final int compareHalf(MutableBigInteger b) {
    int blen = b.intLen;
    int len = intLen;
    if (len <= 0) {
      return blen <= 0 ? 0 : -1;
    }
    if (len > blen) {
      return 1;
    }
    if (len < blen - 1) {
      return -1;
    }
    int[] bval = b.value;
    int bstart = 0;
    int carry = 0;
    // Only 2 cases left:len == blen or len == blen - 1
    if (len != blen) { // len == blen - 1
      if (bval[bstart] == 1) {
        ++bstart;
        carry = 0x80000000;
      } else {
        return -1;
      }
    }
    // compare values with right-shifted values of b,
    // carrying shifted-out bits across words
    int[] val = value;
    for (int i = offset, j = bstart; i < len + offset; ) {
      int bv = bval[j++];
      long hb = ((bv >>> 1) + carry) & LONG_MASK;
      long v = val[i++] & LONG_MASK;
      if (v != hb) {
        return v < hb ? -1 : 1;
      }
      carry = (bv & 1) << 31; // carray will be either 0x80000000 or 0
    }
    return carry == 0 ? 0 : -1;
  }

  /**
   * Return the index of the lowest set bit in this MutableBigInteger. If the
   * magnitude of this MutableBigInteger is zero, -1 is returned.
   */
  private final int getLowestSetBit() {
    if (intLen == 0) {
      return -1;
    }
    int j, b;
    for (j = intLen - 1; (j > 0) && (value[j + offset] == 0); j--) {
      ;
    }
    b = value[j + offset];
    if (b == 0) {
      return -1;
    }
    return ((intLen - 1 - j) << 5) + Integer.numberOfTrailingZeros(b);
  }

  /**
   * Return the int in use in this MutableBigInteger at the specified
   * index. This method is not used because it is not inlined on all
   * platforms.
   */
  private final int getInt(int index) {
    return value[offset + index];
  }

  /**
   * Return a long which is equal to the unsigned value of the int in
   * use in this MutableBigInteger at the specified index. This method is
   * not used because it is not inlined on all platforms.
   */
  private final long getLong(int index) {
    return value[offset + index] & LONG_MASK;
  }

  /**
   * Ensure that the MutableBigInteger is in normal form, specifically
   * making sure that there are no leading zeros, and that if the
   * magnitude is zero, then intLen is zero.
   */
  final void normalize() {
    if (intLen == 0) {
      offset = 0;
      return;
    }

    int index = offset;
    if (value[index] != 0) {
      return;
    }

    int indexBound = index + intLen;
    do {
      index++;
    } while (index < indexBound && value[index] == 0);

    int numZeros = index - offset;
    intLen -= numZeros;
    offset = (intLen == 0 ? 0 : offset + numZeros);
  }

  /**
   * If this MutableBigInteger cannot hold len words, increase the size
   * of the value array to len words.
   */
  private final void ensureCapacity(int len) {
    if (value.length < len) {
      value = new int[len];
      offset = 0;
      intLen = len;
    }
  }

  /**
   * Convert this MutableBigInteger into an int array with no leading
   * zeros, of a length that is equal to this MutableBigInteger's intLen.
   */
  int[] toIntArray() {
    int[] result = new int[intLen];
    for (int i = 0; i < intLen; i++) {
      result[i] = value[offset + i];
    }
    return result;
  }

  /**
   * Sets the int at index+offset in this MutableBigInteger to val.
   * This does not get inlined on all platforms so it is not used
   * as often as originally intended.
   */
  void setInt(int index, int val) {
    value[offset + index] = val;
  }

  /**
   * Sets this MutableBigInteger's value array to the specified array.
   * The intLen is set to the specified length.
   */
  void setValue(int[] val, int length) {
    value = val;
    intLen = length;
    offset = 0;
  }

  /**
   * Sets this MutableBigInteger's value array to a copy of the specified
   * array. The intLen is set to the length of the new array.
   */
  void copyValue(MutableBigInteger src) {
    int len = src.intLen;
    if (value.length < len) {
      value = new int[len];
    }
    System.arraycopy(src.value, src.offset, value, 0, len);
    intLen = len;
    offset = 0;
  }

  /**
   * Sets this MutableBigInteger's value array to a copy of the specified
   * array. The intLen is set to the length of the specified array.
   */
  void copyValue(int[] val) {
    int len = val.length;
    if (value.length < len) {
      value = new int[len];
    }
    System.arraycopy(val, 0, value, 0, len);
    intLen = len;
    offset = 0;
  }

  /**
   * Returns true iff this MutableBigInteger has a value of one.
   */
  boolean isOne() {
    return (intLen == 1) && (value[offset] == 1);
  }

  /**
   * Returns true iff this MutableBigInteger has a value of zero.
   */
  boolean isZero() {
    return (intLen == 0);
  }

  /**
   * Returns true iff this MutableBigInteger is even.
   */
  boolean isEven() {
    return (intLen == 0) || ((value[offset + intLen - 1] & 1) == 0);
  }

  /**
   * Returns true iff this MutableBigInteger is odd.
   */
  boolean isOdd() {
    return isZero() ? false : ((value[offset + intLen - 1] & 1) == 1);
  }

  /**
   * Returns true iff this MutableBigInteger is in normal form. A
   * MutableBigInteger is in normal form if it has no leading zeros
   * after the offset, and intLen + offset <= value.length.
   */
  boolean isNormal() {
    if (intLen + offset > value.length) {
      return false;
    }
    if (intLen == 0) {
      return true;
    }
    return (value[offset] != 0);
  }

  /**
   * Returns a String representation of this MutableBigInteger in radix 10.
   */
  public String toString() {
    BigInteger b = toBigInteger(1);
    return b.toString();
  }

  /**
   * Like {@link #rightShift(int)} but {@code n} can be greater than the length of the number.
   */
  void safeRightShift(int n) {
    if (n / 32 >= intLen) {
      reset();
    } else {
      rightShift(n);
    }
  }

  /**
   * Right shift this MutableBigInteger n bits. The MutableBigInteger is left
   * in normal form.
   */
  void rightShift(int n) {
    if (intLen == 0) {
      return;
    }
    int nInts = n >>> 5;
    int nBits = n & 0x1F;
    this.intLen -= nInts;
    if (nBits == 0) {
      return;
    }
    int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
    if (nBits >= bitsInHighWord) {
      this.primitiveLeftShift(32 - nBits);
      this.intLen--;
    } else {
      primitiveRightShift(nBits);
    }
  }

  /**
   * Like {@link #leftShift(int)} but {@code n} can be zero.
   */
  void safeLeftShift(int n) {
    if (n > 0) {
      leftShift(n);
    }
  }

  /**
   * Left shift this MutableBigInteger n bits.
   */
  void leftShift(int n) {
        /*
         * If there is enough storage space in this MutableBigInteger already
         * the available space will be used. Space to the right of the used
         * ints in the value array is faster to utilize, so the extra space
         * will be taken from the right if possible.
         */
    if (intLen == 0) {
      return;
    }
    int nInts = n >>> 5;
    int nBits = n & 0x1F;
    int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);

    // If shift can be done without moving words, do so
    if (n <= (32 - bitsInHighWord)) {
      primitiveLeftShift(nBits);
      return;
    }

    int newLen = intLen + nInts + 1;
    if (nBits <= (32 - bitsInHighWord)) {
      newLen--;
    }
    if (value.length < newLen) {
      // The array must grow
      int[] result = new int[newLen];
      for (int i = 0; i < intLen; i++) {
        result[i] = value[offset + i];
      }
      setValue(result, newLen);
    } else if (value.length - offset >= newLen) {
      // Use space on right
      for (int i = 0; i < newLen - intLen; i++) {
        value[offset + intLen + i] = 0;
      }
    } else {
      // Must use space on left
      for (int i = 0; i < intLen; i++) {
        value[i] = value[offset + i];
      }
      for (int i = intLen; i < newLen; i++) {
        value[i] = 0;
      }
      offset = 0;
    }
    intLen = newLen;
    if (nBits == 0) {
      return;
    }
    if (nBits <= (32 - bitsInHighWord)) {
      primitiveLeftShift(nBits);
    } else {
      primitiveRightShift(32 - nBits);
    }
  }

  /**
   * A primitive used for division. This method adds in one multiple of the
   * divisor a back to the dividend result at a specified offset. It is used
   * when qhat was estimated too large, and must be adjusted.
   */
  private int divadd(int[] a, int[] result, int offset) {
    long carry = 0;

    for (int j = a.length - 1; j >= 0; j--) {
      long sum = (a[j] & LONG_MASK) +
          (result[j + offset] & LONG_MASK) + carry;
      result[j + offset] = (int) sum;
      carry = sum >>> 32;
    }
    return (int) carry;
  }

  /**
   * This method is used for division. It multiplies an n word input a by one
   * word input x, and subtracts the n word product from q. This is needed
   * when subtracting qhat*divisor from dividend.
   */
  private int mulsub(int[] q, int[] a, int x, int len, int offset) {
    long xLong = x & LONG_MASK;
    long carry = 0;
    offset += len;

    for (int j = len - 1; j >= 0; j--) {
      long product = (a[j] & LONG_MASK) * xLong + carry;
      long difference = q[offset] - product;
      q[offset--] = (int) difference;
      carry = (product >>> 32)
          + (((difference & LONG_MASK) >
          (((~(int) product) & LONG_MASK))) ? 1 : 0);
    }
    return (int) carry;
  }

  /**
   * The method is the same as mulsun, except the fact that q array is not
   * updated, the only result of the method is borrow flag.
   */
  private int mulsubBorrow(int[] q, int[] a, int x, int len, int offset) {
    long xLong = x & LONG_MASK;
    long carry = 0;
    offset += len;
    for (int j = len - 1; j >= 0; j--) {
      long product = (a[j] & LONG_MASK) * xLong + carry;
      long difference = q[offset--] - product;
      carry = (product >>> 32)
          + (((difference & LONG_MASK) >
          (((~(int) product) & LONG_MASK))) ? 1 : 0);
    }
    return (int) carry;
  }

  /**
   * Right shift this MutableBigInteger n bits, where n is
   * less than 32.
   * Assumes that intLen > 0, n > 0 for speed
   */
  private final void primitiveRightShift(int n) {
    int[] val = value;
    int n2 = 32 - n;
    for (int i = offset + intLen - 1, c = val[i]; i > offset; i--) {
      int b = c;
      c = val[i - 1];
      val[i] = (c << n2) | (b >>> n);
    }
    val[offset] >>>= n;
  }

  /**
   * Left shift this MutableBigInteger n bits, where n is
   * less than 32.
   * Assumes that intLen > 0, n > 0 for speed
   */
  private final void primitiveLeftShift(int n) {
    int[] val = value;
    int n2 = 32 - n;
    for (int i = offset, c = val[i], m = i + intLen - 1; i < m; i++) {
      int b = c;
      c = val[i + 1];
      val[i] = (b << n) | (c >>> n2);
    }
    val[offset + intLen - 1] <<= n;
  }

  /**
   * Returns a {@code BigInteger} equal to the {@code n}
   * low ints of this number.
   */
  private BigInteger getLower(int n) {
    if (isZero()) {
      return BigInteger.ZERO;
    } else if (intLen < n) {
      return toBigInteger(1);
    } else {
      // strip zeros
      int len = n;
      while (len > 0 && value[offset + intLen - len] == 0) {
        len--;
      }
      int sign = len > 0 ? 1 : 0;
      return new BigInteger(Arrays.copyOfRange(value, offset + intLen - len, offset + intLen),
          sign);
    }
  }

  /**
   * Discards all ints whose index is greater than {@code n}.
   */
  private void keepLower(int n) {
    if (intLen >= n) {
      offset += intLen - n;
      intLen = n;
    }
  }

  /**
   * Adds the contents of two MutableBigInteger objects.The result
   * is placed within this MutableBigInteger.
   * The contents of the addend are not changed.
   */
  void add(MutableBigInteger addend) {
    int x = intLen;
    int y = addend.intLen;
    int resultLen = (intLen > addend.intLen ? intLen : addend.intLen);
    int[] result = (value.length < resultLen ? new int[resultLen] : value);

    int rstart = result.length - 1;
    long sum;
    long carry = 0;

    // Add common parts of both numbers
    while (x > 0 && y > 0) {
      x--;
      y--;
      sum = (value[x + offset] & LONG_MASK) +
          (addend.value[y + addend.offset] & LONG_MASK) + carry;
      result[rstart--] = (int) sum;
      carry = sum >>> 32;
    }

    // Add remainder of the longer number
    while (x > 0) {
      x--;
      if (carry == 0 && result == value && rstart == (x + offset)) {
        return;
      }
      sum = (value[x + offset] & LONG_MASK) + carry;
      result[rstart--] = (int) sum;
      carry = sum >>> 32;
    }
    while (y > 0) {
      y--;
      sum = (addend.value[y + addend.offset] & LONG_MASK) + carry;
      result[rstart--] = (int) sum;
      carry = sum >>> 32;
    }

    if (carry > 0) { // Result must grow in length
      resultLen++;
      if (result.length < resultLen) {
        int temp[] = new int[resultLen];
        // Result one word longer from carry-out; copy low-order
        // bits into new result.
        System.arraycopy(result, 0, temp, 1, result.length);
        temp[0] = 1;
        result = temp;
      } else {
        result[rstart--] = 1;
      }
    }

    value = result;
    intLen = resultLen;
    offset = result.length - resultLen;
  }

  /**
   * Adds the value of {@code addend} shifted {@code n} ints to the left.
   * Has the same effect as {@code addend.leftShift(32*ints); add(addend);}
   * but doesn't change the value of {@code addend}.
   */
  void addShifted(MutableBigInteger addend, int n) {
    if (addend.isZero()) {
      return;
    }

    int x = intLen;
    int y = addend.intLen + n;
    int resultLen = (intLen > y ? intLen : y);
    int[] result = (value.length < resultLen ? new int[resultLen] : value);

    int rstart = result.length - 1;
    long sum;
    long carry = 0;

    // Add common parts of both numbers
    while (x > 0 && y > 0) {
      x--;
      y--;
      int bval = y + addend.offset < addend.value.length ? addend.value[y + addend.offset] : 0;
      sum = (value[x + offset] & LONG_MASK) +
          (bval & LONG_MASK) + carry;
      result[rstart--] = (int) sum;
      carry = sum >>> 32;
    }

    // Add remainder of the longer number
    while (x > 0) {
      x--;
      if (carry == 0 && result == value && rstart == (x + offset)) {
        return;
      }
      sum = (value[x + offset] & LONG_MASK) + carry;
      result[rstart--] = (int) sum;
      carry = sum >>> 32;
    }
    while (y > 0) {
      y--;
      int bval = y + addend.offset < addend.value.length ? addend.value[y + addend.offset] : 0;
      sum = (bval & LONG_MASK) + carry;
      result[rstart--] = (int) sum;
      carry = sum >>> 32;
    }

    if (carry > 0) { // Result must grow in length
      resultLen++;
      if (result.length < resultLen) {
        int temp[] = new int[resultLen];
        // Result one word longer from carry-out; copy low-order
        // bits into new result.
        System.arraycopy(result, 0, temp, 1, result.length);
        temp[0] = 1;
        result = temp;
      } else {
        result[rstart--] = 1;
      }
    }

    value = result;
    intLen = resultLen;
    offset = result.length - resultLen;
  }

  /**
   * Like {@link #addShifted(MutableBigInteger, int)} but {@code this.intLen} must
   * not be greater than {@code n}. In other words, concatenates {@code this}
   * and {@code addend}.
   */
  void addDisjoint(MutableBigInteger addend, int n) {
    if (addend.isZero()) {
      return;
    }

    int x = intLen;
    int y = addend.intLen + n;
    int resultLen = (intLen > y ? intLen : y);
    int[] result;
    if (value.length < resultLen) {
      result = new int[resultLen];
    } else {
      result = value;
      Arrays.fill(value, offset + intLen, value.length, 0);
    }

    int rstart = result.length - 1;

    // copy from this if needed
    System.arraycopy(value, offset, result, rstart + 1 - x, x);
    y -= x;
    rstart -= x;

    int len = Math.min(y, addend.value.length - addend.offset);
    System.arraycopy(addend.value, addend.offset, result, rstart + 1 - y, len);

    // zero the gap
    for (int i = rstart + 1 - y + len; i < rstart + 1; i++) {
      result[i] = 0;
    }

    value = result;
    intLen = resultLen;
    offset = result.length - resultLen;
  }

  /**
   * Adds the low {@code n} ints of {@code addend}.
   */
  void addLower(MutableBigInteger addend, int n) {
    MutableBigInteger a = new MutableBigInteger(addend);
    if (a.offset + a.intLen >= n) {
      a.offset = a.offset + a.intLen - n;
      a.intLen = n;
    }
    a.normalize();
    add(a);
  }

  /**
   * Subtracts the smaller of this and b from the larger and places the
   * result into this MutableBigInteger.
   */
  int subtract(MutableBigInteger b) {
    MutableBigInteger a = this;

    int[] result = value;
    int sign = a.compare(b);

    if (sign == 0) {
      reset();
      return 0;
    }
    if (sign < 0) {
      MutableBigInteger tmp = a;
      a = b;
      b = tmp;
    }

    int resultLen = a.intLen;
    if (result.length < resultLen) {
      result = new int[resultLen];
    }

    long diff = 0;
    int x = a.intLen;
    int y = b.intLen;
    int rstart = result.length - 1;

    // Subtract common parts of both numbers
    while (y > 0) {
      x--;
      y--;

      diff = (a.value[x + a.offset] & LONG_MASK) -
          (b.value[y + b.offset] & LONG_MASK) - ((int) -(diff >> 32));
      result[rstart--] = (int) diff;
    }
    // Subtract remainder of longer number
    while (x > 0) {
      x--;
      diff = (a.value[x + a.offset] & LONG_MASK) - ((int) -(diff >> 32));
      result[rstart--] = (int) diff;
    }

    value = result;
    intLen = resultLen;
    offset = value.length - resultLen;
    normalize();
    return sign;
  }

  /**
   * Subtracts the smaller of a and b from the larger and places the result
   * into the larger. Returns 1 if the answer is in a, -1 if in b, 0 if no
   * operation was performed.
   */
  private int difference(MutableBigInteger b) {
    MutableBigInteger a = this;
    int sign = a.compare(b);
    if (sign == 0) {
      return 0;
    }
    if (sign < 0) {
      MutableBigInteger tmp = a;
      a = b;
      b = tmp;
    }

    long diff = 0;
    int x = a.intLen;
    int y = b.intLen;

    // Subtract common parts of both numbers
    while (y > 0) {
      x--;
      y--;
      diff = (a.value[a.offset + x] & LONG_MASK) -
          (b.value[b.offset + y] & LONG_MASK) - ((int) -(diff >> 32));
      a.value[a.offset + x] = (int) diff;
    }
    // Subtract remainder of longer number
    while (x > 0) {
      x--;
      diff = (a.value[a.offset + x] & LONG_MASK) - ((int) -(diff >> 32));
      a.value[a.offset + x] = (int) diff;
    }

    a.normalize();
    return sign;
  }

  /**
   * Multiply the contents of two MutableBigInteger objects. The result is
   * placed into MutableBigInteger z. The contents of y are not changed.
   */
  void multiply(MutableBigInteger y, MutableBigInteger z) {
    int xLen = intLen;
    int yLen = y.intLen;
    int newLen = xLen + yLen;

    // Put z into an appropriate state to receive product
    if (z.value.length < newLen) {
      z.value = new int[newLen];
    }
    z.offset = 0;
    z.intLen = newLen;

    // The first iteration is hoisted out of the loop to avoid extra add
    long carry = 0;
    for (int j = yLen - 1, k = yLen + xLen - 1; j >= 0; j--, k--) {
      long product = (y.value[j + y.offset] & LONG_MASK) *
          (value[xLen - 1 + offset] & LONG_MASK) + carry;
      z.value[k] = (int) product;
      carry = product >>> 32;
    }
    z.value[xLen - 1] = (int) carry;

    // Perform the multiplication word by word
    for (int i = xLen - 2; i >= 0; i--) {
      carry = 0;
      for (int j = yLen - 1, k = yLen + i; j >= 0; j--, k--) {
        long product = (y.value[j + y.offset] & LONG_MASK) *
            (value[i + offset] & LONG_MASK) +
            (z.value[k] & LONG_MASK) + carry;
        z.value[k] = (int) product;
        carry = product >>> 32;
      }
      z.value[i] = (int) carry;
    }

    // Remove leading zeros from product
    z.normalize();
  }

  /**
   * Multiply the contents of this MutableBigInteger by the word y. The
   * result is placed into z.
   */
  void mul(int y, MutableBigInteger z) {
    if (y == 1) {
      z.copyValue(this);
      return;
    }

    if (y == 0) {
      z.clear();
      return;
    }

    // Perform the multiplication word by word
    long ylong = y & LONG_MASK;
    int[] zval = (z.value.length < intLen + 1 ? new int[intLen + 1]
        : z.value);
    long carry = 0;
    for (int i = intLen - 1; i >= 0; i--) {
      long product = ylong * (value[i + offset] & LONG_MASK) + carry;
      zval[i + 1] = (int) product;
      carry = product >>> 32;
    }

    if (carry == 0) {
      z.offset = 1;
      z.intLen = intLen;
    } else {
      z.offset = 0;
      z.intLen = intLen + 1;
      zval[0] = (int) carry;
    }
    z.value = zval;
  }

  /**
   * This method is used for division of an n word dividend by a one word
   * divisor. The quotient is placed into quotient. The one word divisor is
   * specified by divisor.
   *
   * @return the remainder of the division is returned.
   *
   */
  int divideOneWord(int divisor, MutableBigInteger quotient) {
    long divisorLong = divisor & LONG_MASK;

    // Special case of one word dividend
    if (intLen == 1) {
      long dividendValue = value[offset] & LONG_MASK;
      int q = (int) (dividendValue / divisorLong);
      int r = (int) (dividendValue - q * divisorLong);
      quotient.value[0] = q;
      quotient.intLen = (q == 0) ? 0 : 1;
      quotient.offset = 0;
      return r;
    }

    if (quotient.value.length < intLen) {
      quotient.value = new int[intLen];
    }
    quotient.offset = 0;
    quotient.intLen = intLen;

    // Normalize the divisor
    int shift = Integer.numberOfLeadingZeros(divisor);

    int rem = value[offset];
    long remLong = rem & LONG_MASK;
    if (remLong < divisorLong) {
      quotient.value[0] = 0;
    } else {
      quotient.value[0] = (int) (remLong / divisorLong);
      rem = (int) (remLong - (quotient.value[0] * divisorLong));
      remLong = rem & LONG_MASK;
    }
    int xlen = intLen;
    while (--xlen > 0) {
      long dividendEstimate = (remLong << 32) |
          (value[offset + intLen - xlen] & LONG_MASK);
      int q;
      if (dividendEstimate >= 0) {
        q = (int) (dividendEstimate / divisorLong);
        rem = (int) (dividendEstimate - q * divisorLong);
      } else {
        long tmp = divWord(dividendEstimate, divisor);
        q = (int) (tmp & LONG_MASK);
        rem = (int) (tmp >>> 32);
      }
      quotient.value[intLen - xlen] = q;
      remLong = rem & LONG_MASK;
    }

    quotient.normalize();
    // Unnormalize
    if (shift > 0) {
      return rem % divisor;
    } else {
      return rem;
    }
  }

  /**
   * Calculates the quotient of this div b and places the quotient in the
   * provided MutableBigInteger objects and the remainder object is returned.
   *
   */
  MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient) {
    return divide(b, quotient, true);
  }

  MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient, boolean needRemainder) {
    if (b.intLen < BigInteger.BURNIKEL_ZIEGLER_THRESHOLD ||
        intLen - b.intLen < BigInteger.BURNIKEL_ZIEGLER_OFFSET) {
      return divideKnuth(b, quotient, needRemainder);
    } else {
      return divideAndRemainderBurnikelZiegler(b, quotient);
    }
  }

  /**
   * @see #divideKnuth(MutableBigInteger, MutableBigInteger, boolean)
   */
  MutableBigInteger divideKnuth(MutableBigInteger b, MutableBigInteger quotient) {
    return divideKnuth(b, quotient, true);
  }

  /**
   * Calculates the quotient of this div b and places the quotient in the
   * provided MutableBigInteger objects and the remainder object is returned.
   *
   * Uses Algorithm D in Knuth section 4.3.1.
   * Many optimizations to that algorithm have been adapted from the Colin
   * Plumb C library.
   * It special cases one word divisors for speed. The content of b is not
   * changed.
   *
   */
  MutableBigInteger divideKnuth(MutableBigInteger b, MutableBigInteger quotient,
      boolean needRemainder) {
    if (b.intLen == 0) {
      throw new ArithmeticException("BigInteger divide by zero");
    }

    // Dividend is zero
    if (intLen == 0) {
      quotient.intLen = quotient.offset = 0;
      return needRemainder ? new MutableBigInteger() : null;
    }

    int cmp = compare(b);
    // Dividend less than divisor
    if (cmp < 0) {
      quotient.intLen = quotient.offset = 0;
      return needRemainder ? new MutableBigInteger(this) : null;
    }
    // Dividend equal to divisor
    if (cmp == 0) {
      quotient.value[0] = quotient.intLen = 1;
      quotient.offset = 0;
      return needRemainder ? new MutableBigInteger() : null;
    }

    quotient.clear();
    // Special case one word divisor
    if (b.intLen == 1) {
      int r = divideOneWord(b.value[b.offset], quotient);
      if (needRemainder) {
        if (r == 0) {
          return new MutableBigInteger();
        }
        return new MutableBigInteger(r);
      } else {
        return null;
      }
    }

    // Cancel common powers of two if we're above the KNUTH_POW2_* thresholds
    if (intLen >= KNUTH_POW2_THRESH_LEN) {
      int trailingZeroBits = Math.min(getLowestSetBit(), b.getLowestSetBit());
      if (trailingZeroBits >= KNUTH_POW2_THRESH_ZEROS * 32) {
        MutableBigInteger a = new MutableBigInteger(this);
        b = new MutableBigInteger(b);
        a.rightShift(trailingZeroBits);
        b.rightShift(trailingZeroBits);
        MutableBigInteger r = a.divideKnuth(b, quotient);
        r.leftShift(trailingZeroBits);
        return r;
      }
    }

    return divideMagnitude(b, quotient, needRemainder);
  }

  /**
   * Computes {@code this/b} and {@code this%b} using the
   * <a href="http://cr.yp.to/bib/1998/burnikel.ps"> Burnikel-Ziegler algorithm</a>.
   * This method implements algorithm 3 from pg. 9 of the Burnikel-Ziegler paper.
   * The parameter beta was chosen to b 2<sup>32</sup> so almost all shifts are
   * multiples of 32 bits.<br/>
   * {@code this} and {@code b} must be nonnegative.
   * @param b the divisor
   * @param quotient output parameter for {@code this/b}
   * @return the remainder
   */
  MutableBigInteger divideAndRemainderBurnikelZiegler(MutableBigInteger b,
      MutableBigInteger quotient) {
    int r = intLen;
    int s = b.intLen;

    // Clear the quotient
    quotient.offset = quotient.intLen = 0;

    if (r < s) {
      return this;
    } else {
      // Unlike Knuth division, we don't check for common powers of two here because
      // BZ already runs faster if both numbers contain powers of two and cancelling them has no
      // additional benefit.

      // step 1: let m = min{2^k | (2^k)*BURNIKEL_ZIEGLER_THRESHOLD > s}
      int m = 1 << (32 - Integer.numberOfLeadingZeros(s / BigInteger.BURNIKEL_ZIEGLER_THRESHOLD));

      int j = (s + m - 1) / m;      // step 2a: j = ceil(s/m)
      int n = j * m;            // step 2b: block length in 32-bit units
      long n32 = 32L * n;         // block length in bits
      int sigma = (int) Math
          .max(0, n32 - b.bitLength());   // step 3: sigma = max{T | (2^T)*B < beta^n}
      MutableBigInteger bShifted = new MutableBigInteger(b);
      bShifted.safeLeftShift(sigma);   // step 4a: shift b so its length is a multiple of n
      MutableBigInteger aShifted = new MutableBigInteger(this);
      aShifted.safeLeftShift(sigma);     // step 4b: shift a by the same amount

      // step 5: t is the number of blocks needed to accommodate a plus one additional bit
      int t = (int) ((aShifted.bitLength() + n32) / n32);
      if (t < 2) {
        t = 2;
      }

      // step 6: conceptually split a into blocks a[t-1], ..., a[0]
      MutableBigInteger a1 = aShifted.getBlock(t - 1, t, n);   // the most significant block of a

      // step 7: z[t-2] = [a[t-1], a[t-2]]
      MutableBigInteger z = aShifted
          .getBlock(t - 2, t, n);    // the second to most significant block
      z.addDisjoint(a1, n);   // z[t-2]

      // do schoolbook division on blocks, dividing 2-block numbers by 1-block numbers
      MutableBigInteger qi = new MutableBigInteger();
      MutableBigInteger ri;
      for (int i = t - 2; i > 0; i--) {
        // step 8a: compute (qi,ri) such that z=b*qi+ri
        ri = z.divide2n1n(bShifted, qi);

        // step 8b: z = [ri, a[i-1]]
        z = aShifted.getBlock(i - 1, t, n);   // a[i-1]
        z.addDisjoint(ri, n);
        quotient.addShifted(qi, i * n);   // update q (part of step 9)
      }
      // final iteration of step 8: do the loop one more time for i=0 but leave z unchanged
      ri = z.divide2n1n(bShifted, qi);
      quotient.add(qi);

      ri.rightShift(sigma);   // step 9: a and b were shifted, so shift back
      return ri;
    }
  }

  /**
   * This method implements algorithm 1 from pg. 4 of the Burnikel-Ziegler paper.
   * It divides a 2n-digit number by a n-digit number.<br/>
   * The parameter beta is 2<sup>32</sup> so all shifts are multiples of 32 bits.
   * <br/>
   * {@code this} must be a nonnegative number such that {@code this.bitLength() <= 2*b.bitLength()}
   * @param b a positive number such that {@code b.bitLength()} is even
   * @param quotient output parameter for {@code this/b}
   * @return {@code this%b}
   */
  private MutableBigInteger divide2n1n(MutableBigInteger b, MutableBigInteger quotient) {
    int n = b.intLen;

    // step 1: base case
    if (n % 2 != 0 || n < BigInteger.BURNIKEL_ZIEGLER_THRESHOLD) {
      return divideKnuth(b, quotient);
    }

    // step 2: view this as [a1,a2,a3,a4] where each ai is n/2 ints or less
    MutableBigInteger aUpper = new MutableBigInteger(this);
    aUpper.safeRightShift(32 * (n / 2));   // aUpper = [a1,a2,a3]
    keepLower(n / 2);   // this = a4

    // step 3: q1=aUpper/b, r1=aUpper%b
    MutableBigInteger q1 = new MutableBigInteger();
    MutableBigInteger r1 = aUpper.divide3n2n(b, q1);

    // step 4: quotient=[r1,this]/b, r2=[r1,this]%b
    addDisjoint(r1, n / 2);   // this = [r1,this]
    MutableBigInteger r2 = divide3n2n(b, quotient);

    // step 5: let quotient=[q1,quotient] and return r2
    quotient.addDisjoint(q1, n / 2);
    return r2;
  }

  /**
   * This method implements algorithm 2 from pg. 5 of the Burnikel-Ziegler paper.
   * It divides a 3n-digit number by a 2n-digit number.<br/>
   * The parameter beta is 2<sup>32</sup> so all shifts are multiples of 32 bits.<br/>
   * <br/>
   * {@code this} must be a nonnegative number such that {@code 2*this.bitLength() <= 3*b.bitLength()}
   * @param quotient output parameter for {@code this/b}
   * @return {@code this%b}
   */
  private MutableBigInteger divide3n2n(MutableBigInteger b, MutableBigInteger quotient) {
    int n = b.intLen / 2;   // half the length of b in ints

    // step 1: view this as [a1,a2,a3] where each ai is n ints or less; let a12=[a1,a2]
    MutableBigInteger a12 = new MutableBigInteger(this);
    a12.safeRightShift(32 * n);

    // step 2: view b as [b1,b2] where each bi is n ints or less
    MutableBigInteger b1 = new MutableBigInteger(b);
    b1.safeRightShift(n * 32);
    BigInteger b2 = b.getLower(n);

    MutableBigInteger r;
    MutableBigInteger d;
    if (compareShifted(b, n) < 0) {
      // step 3a: if a1<b1, let quotient=a12/b1 and r=a12%b1
      r = a12.divide2n1n(b1, quotient);

      // step 4: d=quotient*b2
      d = new MutableBigInteger(quotient.toBigInteger().multiply(b2));
    } else {
      // step 3b: if a1>=b1, let quotient=beta^n-1 and r=a12-b1*2^n+b1
      quotient.ones(n);
      a12.add(b1);
      b1.leftShift(32 * n);
      a12.subtract(b1);
      r = a12;

      // step 4: d=quotient*b2=(b2 << 32*n) - b2
      d = new MutableBigInteger(b2);
      d.leftShift(32 * n);
      d.subtract(new MutableBigInteger(b2));
    }

    // step 5: r = r*beta^n + a3 - d (paper says a4)
    // However, don't subtract d until after the while loop so r doesn't become negative
    r.leftShift(32 * n);
    r.addLower(this, n);

    // step 6: add b until r>=d
    while (r.compare(d) < 0) {
      r.add(b);
      quotient.subtract(MutableBigInteger.ONE);
    }
    r.subtract(d);

    return r;
  }

  /**
   * Returns a {@code MutableBigInteger} containing {@code blockLength} ints from
   * {@code this} number, starting at {@code index*blockLength}.<br/>
   * Used by Burnikel-Ziegler division.
   * @param index the block index
   * @param numBlocks the total number of blocks in {@code this} number
   * @param blockLength length of one block in units of 32 bits
   * @return
   */
  private MutableBigInteger getBlock(int index, int numBlocks, int blockLength) {
    int blockStart = index * blockLength;
    if (blockStart >= intLen) {
      return new MutableBigInteger();
    }

    int blockEnd;
    if (index == numBlocks - 1) {
      blockEnd = intLen;
    } else {
      blockEnd = (index + 1) * blockLength;
    }
    if (blockEnd > intLen) {
      return new MutableBigInteger();
    }

    int[] newVal = Arrays
        .copyOfRange(value, offset + intLen - blockEnd, offset + intLen - blockStart);
    return new MutableBigInteger(newVal);
  }

  /** @see BigInteger#bitLength() */
  long bitLength() {
    if (intLen == 0) {
      return 0;
    }
    return intLen * 32L - Integer.numberOfLeadingZeros(value[offset]);
  }

  /**
   * Internally used  to calculate the quotient of this div v and places the
   * quotient in the provided MutableBigInteger object and the remainder is
   * returned.
   *
   * @return the remainder of the division will be returned.
   */
  long divide(long v, MutableBigInteger quotient) {
    if (v == 0) {
      throw new ArithmeticException("BigInteger divide by zero");
    }

    // Dividend is zero
    if (intLen == 0) {
      quotient.intLen = quotient.offset = 0;
      return 0;
    }
    if (v < 0) {
      v = -v;
    }

    int d = (int) (v >>> 32);
    quotient.clear();
    // Special case on word divisor
    if (d == 0) {
      return divideOneWord((int) v, quotient) & LONG_MASK;
    } else {
      return divideLongMagnitude(v, quotient).toLong();
    }
  }

  private static void copyAndShift(int[] src, int srcFrom, int srcLen, int[] dst, int dstFrom,
      int shift) {
    int n2 = 32 - shift;
    int c = src[srcFrom];
    for (int i = 0; i < srcLen - 1; i++) {
      int b = c;
      c = src[++srcFrom];
      dst[dstFrom + i] = (b << shift) | (c >>> n2);
    }
    dst[dstFrom + srcLen - 1] = c << shift;
  }

  /**
   * Divide this MutableBigInteger by the divisor.
   * The quotient will be placed into the provided quotient object &
   * the remainder object is returned.
   */
  private MutableBigInteger divideMagnitude(MutableBigInteger div,
      MutableBigInteger quotient,
      boolean needRemainder) {
    // assert div.intLen > 1
    // D1 normalize the divisor
    int shift = Integer.numberOfLeadingZeros(div.value[div.offset]);
    // Copy divisor value to protect divisor
    final int dlen = div.intLen;
    int[] divisor;
    MutableBigInteger rem; // Remainder starts as dividend with space for a leading zero
    if (shift > 0) {
      divisor = new int[dlen];
      copyAndShift(div.value, div.offset, dlen, divisor, 0, shift);
      if (Integer.numberOfLeadingZeros(value[offset]) >= shift) {
        int[] remarr = new int[intLen + 1];
        rem = new MutableBigInteger(remarr);
        rem.intLen = intLen;
        rem.offset = 1;
        copyAndShift(value, offset, intLen, remarr, 1, shift);
      } else {
        int[] remarr = new int[intLen + 2];
        rem = new MutableBigInteger(remarr);
        rem.intLen = intLen + 1;
        rem.offset = 1;
        int rFrom = offset;
        int c = 0;
        int n2 = 32 - shift;
        for (int i = 1; i < intLen + 1; i++, rFrom++) {
          int b = c;
          c = value[rFrom];
          remarr[i] = (b << shift) | (c >>> n2);
        }
        remarr[intLen + 1] = c << shift;
      }
    } else {
      divisor = Arrays.copyOfRange(div.value, div.offset, div.offset + div.intLen);
      rem = new MutableBigInteger(new int[intLen + 1]);
      System.arraycopy(value, offset, rem.value, 1, intLen);
      rem.intLen = intLen;
      rem.offset = 1;
    }

    int nlen = rem.intLen;

    // Set the quotient size
    final int limit = nlen - dlen + 1;
    if (quotient.value.length < limit) {
      quotient.value = new int[limit];
      quotient.offset = 0;
    }
    quotient.intLen = limit;
    int[] q = quotient.value;

    // Must insert leading 0 in rem if its length did not change
    if (rem.intLen == nlen) {
      rem.offset = 0;
      rem.value[0] = 0;
      rem.intLen++;
    }

    int dh = divisor[0];
    long dhLong = dh & LONG_MASK;
    int dl = divisor[1];

    // D2 Initialize j
    for (int j = 0; j < limit - 1; j++) {
      // D3 Calculate qhat
      // estimate qhat
      int qhat = 0;
      int qrem = 0;
      boolean skipCorrection = false;
      int nh = rem.value[j + rem.offset];
      int nh2 = nh + 0x80000000;
      int nm = rem.value[j + 1 + rem.offset];

      if (nh == dh) {
        qhat = ~0;
        qrem = nh + nm;
        skipCorrection = qrem + 0x80000000 < nh2;
      } else {
        long nChunk = (((long) nh) << 32) | (nm & LONG_MASK);
        if (nChunk >= 0) {
          qhat = (int) (nChunk / dhLong);
          qrem = (int) (nChunk - (qhat * dhLong));
        } else {
          long tmp = divWord(nChunk, dh);
          qhat = (int) (tmp & LONG_MASK);
          qrem = (int) (tmp >>> 32);
        }
      }

      if (qhat == 0) {
        continue;
      }

      if (!skipCorrection) { // Correct qhat
        long nl = rem.value[j + 2 + rem.offset] & LONG_MASK;
        long rs = ((qrem & LONG_MASK) << 32) | nl;
        long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);

        if (unsignedLongCompare(estProduct, rs)) {
          qhat--;
          qrem = (int) ((qrem & LONG_MASK) + dhLong);
          if ((qrem & LONG_MASK) >= dhLong) {
            estProduct -= (dl & LONG_MASK);
            rs = ((qrem & LONG_MASK) << 32) | nl;
            if (unsignedLongCompare(estProduct, rs)) {
              qhat--;
            }
          }
        }
      }

      // D4 Multiply and subtract
      rem.value[j + rem.offset] = 0;
      int borrow = mulsub(rem.value, divisor, qhat, dlen, j + rem.offset);

      // D5 Test remainder
      if (borrow + 0x80000000 > nh2) {
        // D6 Add back
        divadd(divisor, rem.value, j + 1 + rem.offset);
        qhat--;
      }

      // Store the quotient digit
      q[j] = qhat;
    } // D7 loop on j
    // D3 Calculate qhat
    // estimate qhat
    int qhat = 0;
    int qrem = 0;
    boolean skipCorrection = false;
    int nh = rem.value[limit - 1 + rem.offset];
    int nh2 = nh + 0x80000000;
    int nm = rem.value[limit + rem.offset];

    if (nh == dh) {
      qhat = ~0;
      qrem = nh + nm;
      skipCorrection = qrem + 0x80000000 < nh2;
    } else {
      long nChunk = (((long) nh) << 32) | (nm & LONG_MASK);
      if (nChunk >= 0) {
        qhat = (int) (nChunk / dhLong);
        qrem = (int) (nChunk - (qhat * dhLong));
      } else {
        long tmp = divWord(nChunk, dh);
        qhat = (int) (tmp & LONG_MASK);
        qrem = (int) (tmp >>> 32);
      }
    }
    if (qhat != 0) {
      if (!skipCorrection) { // Correct qhat
        long nl = rem.value[limit + 1 + rem.offset] & LONG_MASK;
        long rs = ((qrem & LONG_MASK) << 32) | nl;
        long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);

        if (unsignedLongCompare(estProduct, rs)) {
          qhat--;
          qrem = (int) ((qrem & LONG_MASK) + dhLong);
          if ((qrem & LONG_MASK) >= dhLong) {
            estProduct -= (dl & LONG_MASK);
            rs = ((qrem & LONG_MASK) << 32) | nl;
            if (unsignedLongCompare(estProduct, rs)) {
              qhat--;
            }
          }
        }
      }

      // D4 Multiply and subtract
      int borrow;
      rem.value[limit - 1 + rem.offset] = 0;
      if (needRemainder) {
        borrow = mulsub(rem.value, divisor, qhat, dlen, limit - 1 + rem.offset);
      } else {
        borrow = mulsubBorrow(rem.value, divisor, qhat, dlen, limit - 1 + rem.offset);
      }

      // D5 Test remainder
      if (borrow + 0x80000000 > nh2) {
        // D6 Add back
        if (needRemainder) {
          divadd(divisor, rem.value, limit - 1 + 1 + rem.offset);
        }
        qhat--;
      }

      // Store the quotient digit
      q[(limit - 1)] = qhat;
    }

    if (needRemainder) {
      // D8 Unnormalize
      if (shift > 0) {
        rem.rightShift(shift);
      }
      rem.normalize();
    }
    quotient.normalize();
    return needRemainder ? rem : null;
  }

  /**
   * Divide this MutableBigInteger by the divisor represented by positive long
   * value. The quotient will be placed into the provided quotient object &
   * the remainder object is returned.
   */
  private MutableBigInteger divideLongMagnitude(long ldivisor, MutableBigInteger quotient) {
    // Remainder starts as dividend with space for a leading zero
    MutableBigInteger rem = new MutableBigInteger(new int[intLen + 1]);
    System.arraycopy(value, offset, rem.value, 1, intLen);
    rem.intLen = intLen;
    rem.offset = 1;

    int nlen = rem.intLen;

    int limit = nlen - 2 + 1;
    if (quotient.value.length < limit) {
      quotient.value = new int[limit];
      quotient.offset = 0;
    }
    quotient.intLen = limit;
    int[] q = quotient.value;

    // D1 normalize the divisor
    int shift = Long.numberOfLeadingZeros(ldivisor);
    if (shift > 0) {
      ldivisor <<= shift;
      rem.leftShift(shift);
    }

    // Must insert leading 0 in rem if its length did not change
    if (rem.intLen == nlen) {
      rem.offset = 0;
      rem.value[0] = 0;
      rem.intLen++;
    }

    int dh = (int) (ldivisor >>> 32);
    long dhLong = dh & LONG_MASK;
    int dl = (int) (ldivisor & LONG_MASK);

    // D2 Initialize j
    for (int j = 0; j < limit; j++) {
      // D3 Calculate qhat
      // estimate qhat
      int qhat = 0;
      int qrem = 0;
      boolean skipCorrection = false;
      int nh = rem.value[j + rem.offset];
      int nh2 = nh + 0x80000000;
      int nm = rem.value[j + 1 + rem.offset];

      if (nh == dh) {
        qhat = ~0;
        qrem = nh + nm;
        skipCorrection = qrem + 0x80000000 < nh2;
      } else {
        long nChunk = (((long) nh) << 32) | (nm & LONG_MASK);
        if (nChunk >= 0) {
          qhat = (int) (nChunk / dhLong);
          qrem = (int) (nChunk - (qhat * dhLong));
        } else {
          long tmp = divWord(nChunk, dh);
          qhat = (int) (tmp & LONG_MASK);
          qrem = (int) (tmp >>> 32);
        }
      }

      if (qhat == 0) {
        continue;
      }

      if (!skipCorrection) { // Correct qhat
        long nl = rem.value[j + 2 + rem.offset] & LONG_MASK;
        long rs = ((qrem & LONG_MASK) << 32) | nl;
        long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);

        if (unsignedLongCompare(estProduct, rs)) {
          qhat--;
          qrem = (int) ((qrem & LONG_MASK) + dhLong);
          if ((qrem & LONG_MASK) >= dhLong) {
            estProduct -= (dl & LONG_MASK);
            rs = ((qrem & LONG_MASK) << 32) | nl;
            if (unsignedLongCompare(estProduct, rs)) {
              qhat--;
            }
          }
        }
      }

      // D4 Multiply and subtract
      rem.value[j + rem.offset] = 0;
      int borrow = mulsubLong(rem.value, dh, dl, qhat, j + rem.offset);

      // D5 Test remainder
      if (borrow + 0x80000000 > nh2) {
        // D6 Add back
        divaddLong(dh, dl, rem.value, j + 1 + rem.offset);
        qhat--;
      }

      // Store the quotient digit
      q[j] = qhat;
    } // D7 loop on j

    // D8 Unnormalize
    if (shift > 0) {
      rem.rightShift(shift);
    }

    quotient.normalize();
    rem.normalize();
    return rem;
  }

  /**
   * A primitive used for division by long.
   * Specialized version of the method divadd.
   * dh is a high part of the divisor, dl is a low part
   */
  private int divaddLong(int dh, int dl, int[] result, int offset) {
    long carry = 0;

    long sum = (dl & LONG_MASK) + (result[1 + offset] & LONG_MASK);
    result[1 + offset] = (int) sum;

    sum = (dh & LONG_MASK) + (result[offset] & LONG_MASK) + carry;
    result[offset] = (int) sum;
    carry = sum >>> 32;
    return (int) carry;
  }

  /**
   * This method is used for division by long.
   * Specialized version of the method sulsub.
   * dh is a high part of the divisor, dl is a low part
   */
  private int mulsubLong(int[] q, int dh, int dl, int x, int offset) {
    long xLong = x & LONG_MASK;
    offset += 2;
    long product = (dl & LONG_MASK) * xLong;
    long difference = q[offset] - product;
    q[offset--] = (int) difference;
    long carry = (product >>> 32)
        + (((difference & LONG_MASK) >
        (((~(int) product) & LONG_MASK))) ? 1 : 0);
    product = (dh & LONG_MASK) * xLong + carry;
    difference = q[offset] - product;
    q[offset--] = (int) difference;
    carry = (product >>> 32)
        + (((difference & LONG_MASK) >
        (((~(int) product) & LONG_MASK))) ? 1 : 0);
    return (int) carry;
  }

  /**
   * Compare two longs as if they were unsigned.
   * Returns true iff one is bigger than two.
   */
  private boolean unsignedLongCompare(long one, long two) {
    return (one + Long.MIN_VALUE) > (two + Long.MIN_VALUE);
  }

  /**
   * This method divides a long quantity by an int to estimate
   * qhat for two multi precision numbers. It is used when
   * the signed value of n is less than zero.
   * Returns long value where high 32 bits contain remainder value and
   * low 32 bits contain quotient value.
   */
  static long divWord(long n, int d) {
    long dLong = d & LONG_MASK;
    long r;
    long q;
    if (dLong == 1) {
      q = (int) n;
      r = 0;
      return (r << 32) | (q & LONG_MASK);
    }

    // Approximate the quotient and remainder
    q = (n >>> 1) / (dLong >>> 1);
    r = n - q * dLong;

    // Correct the approximation
    while (r < 0) {
      r += dLong;
      q--;
    }
    while (r >= dLong) {
      r -= dLong;
      q++;
    }
    // n - q*dlong == r && 0 <= r <dLong, hence we're done.
    return (r << 32) | (q & LONG_MASK);
  }

  /**
   * Calculate GCD of this and b. This and b are changed by the computation.
   */
  MutableBigInteger hybridGCD(MutableBigInteger b) {
    // Use Euclid's algorithm until the numbers are approximately the
    // same length, then use the binary GCD algorithm to find the GCD.
    MutableBigInteger a = this;
    MutableBigInteger q = new MutableBigInteger();

    while (b.intLen != 0) {
      if (Math.abs(a.intLen - b.intLen) < 2) {
        return a.binaryGCD(b);
      }

      MutableBigInteger r = a.divide(b, q);
      a = b;
      b = r;
    }
    return a;
  }

  /**
   * Calculate GCD of this and v.
   * Assumes that this and v are not zero.
   */
  private MutableBigInteger binaryGCD(MutableBigInteger v) {
    // Algorithm B from Knuth section 4.5.2
    MutableBigInteger u = this;
    MutableBigInteger r = new MutableBigInteger();

    // step B1
    int s1 = u.getLowestSetBit();
    int s2 = v.getLowestSetBit();
    int k = (s1 < s2) ? s1 : s2;
    if (k != 0) {
      u.rightShift(k);
      v.rightShift(k);
    }

    // step B2
    boolean uOdd = (k == s1);
    MutableBigInteger t = uOdd ? v : u;
    int tsign = uOdd ? -1 : 1;

    int lb;
    while ((lb = t.getLowestSetBit()) >= 0) {
      // steps B3 and B4
      t.rightShift(lb);
      // step B5
      if (tsign > 0) {
        u = t;
      } else {
        v = t;
      }

      // Special case one word numbers
      if (u.intLen < 2 && v.intLen < 2) {
        int x = u.value[u.offset];
        int y = v.value[v.offset];
        x = binaryGcd(x, y);
        r.value[0] = x;
        r.intLen = 1;
        r.offset = 0;
        if (k > 0) {
          r.leftShift(k);
        }
        return r;
      }

      // step B6
      if ((tsign = u.difference(v)) == 0) {
        break;
      }
      t = (tsign >= 0) ? u : v;
    }

    if (k > 0) {
      u.leftShift(k);
    }
    return u;
  }

  /**
   * Calculate GCD of a and b interpreted as unsigned integers.
   */
  static int binaryGcd(int a, int b) {
    if (b == 0) {
      return a;
    }
    if (a == 0) {
      return b;
    }

    // Right shift a & b till their last bits equal to 1.
    int aZeros = Integer.numberOfTrailingZeros(a);
    int bZeros = Integer.numberOfTrailingZeros(b);
    a >>>= aZeros;
    b >>>= bZeros;

    int t = (aZeros < bZeros ? aZeros : bZeros);

    while (a != b) {
      if ((a + 0x80000000) > (b + 0x80000000)) {  // a > b as unsigned
        a -= b;
        a >>>= Integer.numberOfTrailingZeros(a);
      } else {
        b -= a;
        b >>>= Integer.numberOfTrailingZeros(b);
      }
    }
    return a << t;
  }

  /**
   * Returns the modInverse of this mod p. This and p are not affected by
   * the operation.
   */
  MutableBigInteger mutableModInverse(MutableBigInteger p) {
    // Modulus is odd, use Schroeppel's algorithm
    if (p.isOdd()) {
      return modInverse(p);
    }

    // Base and modulus are even, throw exception
    if (isEven()) {
      throw new ArithmeticException("BigInteger not invertible.");
    }

    // Get even part of modulus expressed as a power of 2
    int powersOf2 = p.getLowestSetBit();

    // Construct odd part of modulus
    MutableBigInteger oddMod = new MutableBigInteger(p);
    oddMod.rightShift(powersOf2);

    if (oddMod.isOne()) {
      return modInverseMP2(powersOf2);
    }

    // Calculate 1/a mod oddMod
    MutableBigInteger oddPart = modInverse(oddMod);

    // Calculate 1/a mod evenMod
    MutableBigInteger evenPart = modInverseMP2(powersOf2);

    // Combine the results using Chinese Remainder Theorem
    MutableBigInteger y1 = modInverseBP2(oddMod, powersOf2);
    MutableBigInteger y2 = oddMod.modInverseMP2(powersOf2);

    MutableBigInteger temp1 = new MutableBigInteger();
    MutableBigInteger temp2 = new MutableBigInteger();
    MutableBigInteger result = new MutableBigInteger();

    oddPart.leftShift(powersOf2);
    oddPart.multiply(y1, result);

    evenPart.multiply(oddMod, temp1);
    temp1.multiply(y2, temp2);

    result.add(temp2);
    return result.divide(p, temp1);
  }

  /*
   * Calculate the multiplicative inverse of this mod 2^k.
   */
  MutableBigInteger modInverseMP2(int k) {
    if (isEven()) {
      throw new ArithmeticException("Non-invertible. (GCD != 1)");
    }

    if (k > 64) {
      return euclidModInverse(k);
    }

    int t = inverseMod32(value[offset + intLen - 1]);

    if (k < 33) {
      t = (k == 32 ? t : t & ((1 << k) - 1));
      return new MutableBigInteger(t);
    }

    long pLong = (value[offset + intLen - 1] & LONG_MASK);
    if (intLen > 1) {
      pLong |= ((long) value[offset + intLen - 2] << 32);
    }
    long tLong = t & LONG_MASK;
    tLong = tLong * (2 - pLong * tLong);  // 1 more Newton iter step
    tLong = (k == 64 ? tLong : tLong & ((1L << k) - 1));

    MutableBigInteger result = new MutableBigInteger(new int[2]);
    result.value[0] = (int) (tLong >>> 32);
    result.value[1] = (int) tLong;
    result.intLen = 2;
    result.normalize();
    return result;
  }

  /**
   * Returns the multiplicative inverse of val mod 2^32.  Assumes val is odd.
   */
  static int inverseMod32(int val) {
    // Newton's iteration!
    int t = val;
    t *= 2 - val * t;
    t *= 2 - val * t;
    t *= 2 - val * t;
    t *= 2 - val * t;
    return t;
  }

  /**
   * Calculate the multiplicative inverse of 2^k mod mod, where mod is odd.
   */
  static MutableBigInteger modInverseBP2(MutableBigInteger mod, int k) {
    // Copy the mod to protect original
    return fixup(new MutableBigInteger(1), new MutableBigInteger(mod), k);
  }

  /**
   * Calculate the multiplicative inverse of this mod mod, where mod is odd.
   * This and mod are not changed by the calculation.
   *
   * This method implements an algorithm due to Richard Schroeppel, that uses
   * the same intermediate representation as Montgomery Reduction
   * ("Montgomery Form").  The algorithm is described in an unpublished
   * manuscript entitled "Fast Modular Reciprocals."
   */
  private MutableBigInteger modInverse(MutableBigInteger mod) {
    MutableBigInteger p = new MutableBigInteger(mod);
    MutableBigInteger f = new MutableBigInteger(this);
    MutableBigInteger g = new MutableBigInteger(p);
    SignedMutableBigInteger c = new SignedMutableBigInteger(1);
    SignedMutableBigInteger d = new SignedMutableBigInteger();
    MutableBigInteger temp = null;
    SignedMutableBigInteger sTemp = null;

    int k = 0;
    // Right shift f k times until odd, left shift d k times
    if (f.isEven()) {
      int trailingZeros = f.getLowestSetBit();
      f.rightShift(trailingZeros);
      d.leftShift(trailingZeros);
      k = trailingZeros;
    }

    // The Almost Inverse Algorithm
    while (!f.isOne()) {
      // If gcd(f, g) != 1, number is not invertible modulo mod
      if (f.isZero()) {
        throw new ArithmeticException("BigInteger not invertible.");
      }

      // If f < g exchange f, g and c, d
      if (f.compare(g) < 0) {
        temp = f;
        f = g;
        g = temp;
        sTemp = d;
        d = c;
        c = sTemp;
      }

      // If f == g (mod 4)
      if (((f.value[f.offset + f.intLen - 1] ^
          g.value[g.offset + g.intLen - 1]) & 3) == 0) {
        f.subtract(g);
        c.signedSubtract(d);
      } else { // If f != g (mod 4)
        f.add(g);
        c.signedAdd(d);
      }

      // Right shift f k times until odd, left shift d k times
      int trailingZeros = f.getLowestSetBit();
      f.rightShift(trailingZeros);
      d.leftShift(trailingZeros);
      k += trailingZeros;
    }

    while (c.sign < 0) {
      c.signedAdd(p);
    }

    return fixup(c, p, k);
  }

  /**
   * The Fixup Algorithm
   * Calculates X such that X = C * 2^(-k) (mod P)
   * Assumes C<P and P is odd.
   */
  static MutableBigInteger fixup(MutableBigInteger c, MutableBigInteger p,
      int k) {
    MutableBigInteger temp = new MutableBigInteger();
    // Set r to the multiplicative inverse of p mod 2^32
    int r = -inverseMod32(p.value[p.offset + p.intLen - 1]);

    for (int i = 0, numWords = k >> 5; i < numWords; i++) {
      // V = R * c (mod 2^j)
      int v = r * c.value[c.offset + c.intLen - 1];
      // c = c + (v * p)
      p.mul(v, temp);
      c.add(temp);
      // c = c / 2^j
      c.intLen--;
    }
    int numBits = k & 0x1f;
    if (numBits != 0) {
      // V = R * c (mod 2^j)
      int v = r * c.value[c.offset + c.intLen - 1];
      v &= ((1 << numBits) - 1);
      // c = c + (v * p)
      p.mul(v, temp);
      c.add(temp);
      // c = c / 2^j
      c.rightShift(numBits);
    }

    // In theory, c may be greater than p at this point (Very rare!)
    while (c.compare(p) >= 0) {
      c.subtract(p);
    }

    return c;
  }

  /**
   * Uses the extended Euclidean algorithm to compute the modInverse of base
   * mod a modulus that is a power of 2. The modulus is 2^k.
   */
  MutableBigInteger euclidModInverse(int k) {
    MutableBigInteger b = new MutableBigInteger(1);
    b.leftShift(k);
    MutableBigInteger mod = new MutableBigInteger(b);

    MutableBigInteger a = new MutableBigInteger(this);
    MutableBigInteger q = new MutableBigInteger();
    MutableBigInteger r = b.divide(a, q);

    MutableBigInteger swapper = b;
    // swap b & r
    b = r;
    r = swapper;

    MutableBigInteger t1 = new MutableBigInteger(q);
    MutableBigInteger t0 = new MutableBigInteger(1);
    MutableBigInteger temp = new MutableBigInteger();

    while (!b.isOne()) {
      r = a.divide(b, q);

      if (r.intLen == 0) {
        throw new ArithmeticException("BigInteger not invertible.");
      }

      swapper = r;
      a = swapper;

      if (q.intLen == 1) {
        t1.mul(q.value[q.offset], temp);
      } else {
        q.multiply(t1, temp);
      }
      swapper = q;
      q = temp;
      temp = swapper;
      t0.add(q);

      if (a.isOne()) {
        return t0;
      }

      r = b.divide(a, q);

      if (r.intLen == 0) {
        throw new ArithmeticException("BigInteger not invertible.");
      }

      swapper = b;
      b = r;

      if (q.intLen == 1) {
        t0.mul(q.value[q.offset], temp);
      } else {
        q.multiply(t0, temp);
      }
      swapper = q;
      q = temp;
      temp = swapper;

      t1.add(q);
    }
    mod.subtract(t1);
    return mod;
  }
}
